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ABSTRACT 

A  three-dimensional  (3D)  coupled  normal  mode  model  for  studying 

sound  propagation  in  a  complex  coastal  environment  is  developed. 
This  development  corresponds  to  a  significant  upgrade  of  an  earlier 
version  of  the  model  in  which  a  flat,  rigid  bottom  was  used.  By 
imposing  the  general  boundary  conditions  for  an  irregular,  non- 
rigid  bottom,  the  coupling  coefficient  integrals  in  the  system  of 
differential  equations  governing  the  mode  amplitude  are  re- 
formulated. The  model  upgrade  entails  a  numerical  implementation  of 
the  revised  formulae.  With  the  improved  physics,  this  latest 
version  is  capable  of  modeling  the  3D  acoustic  wave-field  in 
shallow  water  where  sound  speed,  water  depth  and  sediment 
properties  can  vary  with  horizontal  location.  To  demonstrate  this 
enhanced  capability,  the  model  is  used  here  to  simulate  the 
interactions  of  the  normal  modes  as  they  propagate  up  a  sloping 
bottom. 
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I .     INTRODUCTION 

A .   BACKGROUND 

There  are  three  approaches  to  model  three-dimensional (3D) 
sound  propagation  in  the  ocean:  ray  theory,  parabolic 
equation  approximation  and   normal  mode  theory. 

Ray  theory  gives  an  approximate,  planewave- like  solution 
to  the  wave  equation,  which  is  valid  at  high  enough 
frequencies  and  in  media  with  gradual  variations.  The  ray 
solution  is  constructed  by  raytracing.  The  acoustic  rays 
provide  for  a  visual,  physical  description  of  sound 
transmission  in  the  ocean.  The  ray  solution,  however, 
neglects  sound  diffraction  and  dispersion  and  thus  needs 
corrections  near  caustics  and  turning  points.  These 
corrections  may  sometimes  be  mathematically  complicated.  The 
Hamiltonian  Acoustic  Ray  Tracing  Program  for  the  Ocean  (HARPO) 
is  the  only  3D  ray  theory  model  available  today.  This  computer 
code  was  originally  developed  by  Jones  et  al .  [Ref  .  1]  for  the 
computation  of  3D  rays. 

The  parabolic  equation  approximation  method  (PE)  was 
introduced  by  Tappert  [Ref.  2]  .  PE  is  a  "full -wave"  method 
that  accounts  for  both  sound  diffraction  and  dispersion.  It 
provides  for  numerical  solutions  to  the  wave  equation  which 
are  accurate  for  energy  propagating  at  low  grazing  angles.  The 


accuracy  generally  degrades  as  the  angle  increases.  The 
backscattered  energy  is  generally  neglected  in  this 
approximation.  A  versatile  3D  PE  model  has  been  developed  by 
Lee  et  al .  [Ref.  3]  using  an  implicit  finite  difference 
scheme.  Another  3D  PE  model  was  developed  earlier  by  Baer 
[Ref.  4]  which  uses  a  split-step  Fourier  algorithm.  The  PE 
model  of  Lee  has  a  wider  angle  capability,  i.e.,  it  models 
sound  energy  travelling  at  steeper  angles  more  accurately. 

Finally,  normal  mode  theory  describes  sound  propagation  as 
a  collection  of  eigenf unctions,  called  normal  modes,  which  are 
a  natural  set  of  vertical  vibration  modes  in  the  sound 
channel.  Just  like  PE,  normal  mode  theory  is  a  "full-wave" 
approach.  The  original  normal  mode  theory  assumes  a 
horizontally  stratified  propagation  medium.  This  assumption  is 
valid  for  many  short-range,  deep-water  cases,  where  range  and 
azimuthal  variations  are  negligibly  small.  Pierce  [Ref.  5] 
extended  the  theory  to  account  for  horizontal  sound  speed, 
bathymetry  and  bottom-property  variations.  These  variations 
produce  mode  coupling  phenomena  (in  which  energy  exchange 
between  modes  takes  place) .  A  3D  coupled  normal  mode  model  has 
been  developed  by  Chiu  and  Ehret  [Ref.  6]  .  This  model  is 
capable  of  simulating  mode-mode  interactions  due  to  a  3D 
varying  sound  speed  field.  The  effects  of  bottom  bathymetry 
variations  and  sediment  property,  however,  are  not  modeled. 


B.   THESIS  OBJECTIVES  AND  OUTLINE 

The  main  objective  of  this  thesis  is  to  improve  the  Chiu- 
Ehret  [Ref .  6]  3D  coupled  normal  mode  model  by  including  the 
effects  of  bathymetry  variations  and  sediment  properties  on 
sound  propagation.  The  improved  model  is  useful  for  studying 
sound  propagation  in  shallow  water  environments  where 
significant  bottom  interaction  is  expected. 

In  Chapter  II,  3D  coupled  mode  theory  is  first  reviewed. 
The  formulae  for  the  mode  coupling  coefficients  in  the  system 
of  differential  equations  governing  the  mode  amplitude 
functions  are  derived.  In  the  derivation,  the  general  boundary 
conditions  for  an  irregular,  non-rigid  bottom  are  used. 

In  Chapter  III,  alternative  expressions  for  the  mode 
coupling  coefficients  are  derived.  These  expressions  allow  for 
an  easier  numerical  implementation.  The  improved  model  is  used 
to  examine  the  effects  of  a  sloping  bottom  on  upslope  sound 
propagation.  The  validity  of  the  adiabatic  approximation  is 
also  examined.  Conclusions  are  given  in  Chapter  IV. 

The  products  coming  out  of  this  thesis  are  computer 
subroutines  to  include  bathymetry  variations  and  sediment 
properties  in  the  3D  coupled  mode  model  of  Chiu  and  Ehret 
[Ref.  6] .  The  new  routines  are  listed  in  the  Appendix. 


II.     3D  COUPLED  NORMAL  MODE  THEORY 

In  the  mathematical  formulation  that  follows,  a 
cylindrical  coordinate  system  will  be  used  (see  Fig.  1) .  The 
z-axis  is  perpendicular  to  the  ocean  surface  and  is  positive 
downward,  r  is  range  from  the  source  location  (i.e.,  the 
origin)  and  6  is  the  azimuthal  angle  (positive  clockwise) . 
Sound  speed  in  the  water  column,  c,,  is  a  function  of  r,  z  and 
9,  where  sound  speed  in  the  sediment,  c2,  is  assumed  to  be  r 
and  6  dependent  only.  The  density  of  the  water  column,  p,,  is 
considered  to  be  constant.  The  density  of  the  sediment,  p2,  is 
also  considered  to  be  constant.  The  water- sediment  interface 
is  located  at  z=H(r,0). 

A.   THE  MATHEMATICAL  PROBLEM 

In  the  case  of  isodensity  layers,  the  3D,  homogeneous, 
monof requency  Helmholtz  Equation  governing  the  acoustic 
pressure,  p,  is: 

V2p(z,r,6)  +  k2(z,r,e)p(z,r,e)    =0  (D 

where  k  (z,  r,  0)  =a>/c  (z,  r,  6)  is  the  acoustic  wavenumber,  u  is  the 
source  angular  frequency  and  c  is  sound  speed  (ct  in  the  water 
layer  and  c2  in  the  sediment  layer)  . 

A  quasi -separable  solution  to  Eq.  (1)  is  postulated,  which 
is  locally  a  linear  combination  of  normal  modes  or  depth 
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Figure  1 .    The  model  geometry 
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functions,  Zn 


p(ztr,Q)    =  ^Rn(r,e)Zn(z;r,Q) 


(2) 

n 


where  Rn  are  the  mode  amplitude  functions  and  n  is  the  mode 
number. 

The  normal  modes  Zn  are  required  to  satisfy  the  depth 
equation: 

[-^-  +  k2(z,r,G)     -  k2n(r,e)]Zn(z;r,Q)    =   0        (3) 
oz~ 

where  kn  is  the  horizontal   component  of   the  wavenumber 
(eigenvalue)  associated  with  the  n"1  mode. 

It  can  be  easily  shown,  using  the  boundary  conditions  for 
Zn  (to  be  derived  next)  and  the  depth  equation  (Eq.  (3))  ,  that 
the  normal  modes  form  a  complete  set  of  orthogonal  functions, 
with  the  inverse  of  the  medium  density  p  as  a  weighting 
function  in  the  normalization: 


oo 

\-^Zn^;i,Q)Zm{z;r^)dz  =   5M  (4) 


where  6nm   is  the  Kronecker  delta.  Note  that  the  integration  is 
carried  over  the  entire  depth  from  0  to  oo.  Also,  note  that  the 


density  p  in  the  model  is  considered  to  be  constant  in  each 
layer  (p,  in  the  water  and  p2  in  the  sediment)  . 

B.   THE  BOUNDARY  CONDITIONS 

The  appropriate  boundary  condition  for  the  acoustic 
pressure  at  the  sea  surface  is 

p1(z=0;r,6)  =  0  (5) 

where  the  subscript  1  denotes  the  water  column.  This  pressure 
release  condition  implies  that  the  normal  modes  Zn  must  also 
be  zero  at  the  sea  surface,  i.e., 

Zn(z=0;r,6)  =  0  (6) 

It  also  implies  that  the  horizontal  derivatives  of  Zn  at  z  =  0 
are   zero,  i.e., 

dZn(z=0r,Q)     =  Q  (?) 

dr 


1  azn(z=0;r,9) 

r  ae 


At  the  interface  between  the  sediment  and  the  water 
column,  i.e.,  at  z=H(r,0),  the  boundary  conditions  are 
continuity  of  pressure  and  continuity  of  the  particle  velocity 
component  normal  to  the  interface: 


p1(z=H;r,B)    =  p2(z=H;r,Q)  (9) 


J-VPl.n  =  ±Vp2.h  (10) 

Pi       P2 


The  subscripts  1  and  2  denote  the  water  column  and  the 
sediment,  respectively. 

The  unit  directional  vector  n  normal  to  the  bottom 
interface  is 


A  dH(r,Q)  A  i  6H(r/6)^ 

5}  _   v(z-//(r>e))  ar         r      ae  (id 

IV(z-/f(r,6)l    f  (dH(r,Q)  )2  ,  1  a//(r,6)  J1/2 

ar  r      ae 


where  z,r  and  0  are  the  unit  directional  vectors  associated 
with  the  z,  r  and  6  directions,  respectively.  In  the  case  of 
a  small  bottom  slope,  the  boundary  condition  of  Eq.  (10)  can 
be  approximated  by 


1  dp1(z=H;r,d)  i  dp2(z=H;r,Q) 

px  dz  p2      dz 


The  small  slope  approximation  is  accurate  when 


(12) 


Ite<r'e'l<l  ,  |A3g'f/6>l<l  (13) 

dr  r         36 


Following  Eq.  (9)  and  Eq.  (12),  we  obtain  the  following 
boundary  conditions  at  the  water- sediment  interface  for  the 
normal  modes : 


1  dZln(z=H;r,d)  i  dZ2n(z=H;r,Q] 


p,       dz  p,       dz 


(14) 


Zln(z=H;r,6)    =  Z2n(z=H,r,d)  (15) 

These  boundary  conditions  hold  for  each  individual  normal  mode 
because  they  are  orthogonal  functions. 

The  boundary  condition  for  p  at   z  -*  oo  is 

p2(z-°>,rtQ)    =   0  (16) 

This  implies  that  the  'normal  modes  and  their  horizontal 
derivatives  are  also  zero  as  z  -»  oo: 

Z2n(z^;rfQ)    =0  (17) 


dZ2n(z^;r,Q) 
dr 


=   0  (18) 


1dZ2n(z^rlQ)      ■ 


C.   MODE  COUPLING  COEFFICIENTS  DERIVATION 

Substituting  Eq.  (2)  into  Eq .  (1),  multiplying  by 
Zm(z;r,6)/p,  integrating  over  the  entire  depth  and  finally 
rearranging  terms,  we  obtain  the  coupled  mode  equations 
governing  the  mode  amplitude  functions: 


[Vh*+k\(r,Q)]   RmU,Q)    =  -J2[ABmUie)Vi]Rn(rlQ)+Bmn(rlQ)RnUtQ)] 

n 

(20) 


where  the  two  mode  coupling  coefficients  are  defined  as 

oo 

A^ir,^)    =  2[±Zin(z;r,e)VhZn(z;rlQ)dz  (21) 

o  P 


and 


^(r,6)  =  fj.Zn{z;r,Q)VhZa{z;z,e)dz  (22) 

0 


Vh  is  the  horizontal  gradient  operator,  i.e., 
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V    r  8     A1    d  <23> 


The  mode  coupling  coefficients  are  measures  of  the 
significance  of  exchange  of  acoustic  energy  between  modes 
resulting  from  horizontal  variations  in  the  medium.  As  the 
variations  become  stronger,  the  coupling  coefficients  become 
larger  and  so  is  the  energy  exchange.  In  the  case  of  a 
completely  range  independent  medium,  the  coupling 
coefficients  are  identically  zero  and  the  RHS  of  Eq.  (20) 
vanishes.  In  such  case,  the  modes  propagate  independently  of 
each  other. For  range -dependent  cases,  the  neglect  of  mode 
coupling  leads  to  the  adiabatic  approximation  [Ref .  5] . 

Cylindrical  spreading  can  be  removed  from  the  coupled  mode 
equation  (Eq.  (20)  )  by  replacing  the  mode  amplitude  function 
Rn(r,0)  with  Pn  (r,  8)  /r1/2.  The  result  is 


[  —  +Jc2(r,6)  +  -L  JL  +  _^]p  (r,e)  = 
dr2      m  r2   602  4r2  n 

(24) 

=  Eti/nn(r,e).[VA(r,e)-rP-(r;e)]+Bmn(r,e)Pn(r,e)) 
n  2.r 
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III.   THE  NUMERICAL  MODEL  AND  EXAMPLE  RUNS 

In  this  chapter,  the  procedure  to  upgrade  the  Chiu-Ehret 
model  [Ref.  6]  is  discussed.  The  upgrade  has  entailed  the 
derivation  of  alternative  expressions  for  the  mode  coupling 
coefficients  and  the  generation  of  new  code  to  compute  these 
coefficients  based  on  the  alternative  expressions. 

The  numerical  results  from  two  simple  example  model  runs 
associated  with  two  different  bottom  slopes  are  also  presented 
in  this  chapter.  Both  cases  deal  with  upslope  propagation  in 
isospeed  layers.  These  runs  have  allowed  for  an  examination  of 
the  coupling  between  modes  caused  by  bathymetry  change.  In 
addition,  they  have  allowed  for  an  examination  of  the  validity 
of  the  adiabatic  approximation. 

A.   THE  CHIU-EHRET  APPROACH 

In  the  far  field,  i.e.,  kr>>l  ,  the  coupled  mode  equation 
(Eq.  (24))  for  the  mode  amplitude  functions,  can  be  recast  as 

r>2  1    r)2 

[  — +*-2(r/e)+-i--£-]P_(r/e)  = 

dr2    m  r2  ae2 

(25) 

EGB(r,e),V1Pa(rle)+BB(r/e)Pa(r,e)l 

In  the  Chiu-Ehret  model  [Ref.  6],  Pn   is  decomposed  as 
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Pn(r,Q)    =  Un(r,Q)    ej*»(r'e) 

r  (26) 

4>n(r,6)  =  jkn(r,Q)    dr 


where  Un  is  the  slowly  varying  complex  envelope  of  Pn 
modulating  the  rapidly  varying  two-dimensional  (2D)  adiabatic 
solution,  i.e.,  exp(j0n),  and  #n  is  the  adiabatic  phase.  The 
Chiu-Ehret  model  numerically  computes  the  envelopes  Un  instead 
of  Pn  using  Runge-Kutta  schemes. 

B.   ALTERNATIVE  EXPRESSIONS  FOR  THE  MODE  COUPLING  COEFFICIENTS 

For  simpler  numerical  implementation,  the  expressions  of 
the  mode  coupling  coefficients  in  Eq.  (21)  and  Eq.  (22)  are 
rewritten  in  alternative  forms.  These  alternative  forms  do  not 
require  integrations  of  expressions  involving  the  horizontal 
derivatives  of  normal  modes.  In  the  following,  the  derivation 
of  these  alternative  forms  is  presented. 

1.   VECTOR  MODE  COUPLING  COEFFICIENT,  A^ 

a.  Case  of  nstn 

Applying  the  horizontal  gradient  operator  Vh  to 
both  sides  of  the  depth  equation  (Eq.  (3)),  we  get 
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dz2 


+  2[k(z,r,B)Vhk(z,r,B)    -  kn(r,e)Vhkn(rfO)]Zn   + 


[Jc2(z,r,0)  -  £2n(r,0)]VhZn  =  0 


(27) 


Multiplying  Eq.  (27)  by  Zm(z;r,0)/p  and  then  integrating  over 
the  entire  depth,  we  get 


1  ,   ^A 


rizm £-2  dz  +  f  (k2-kn2)  ±Zn  VhZn  dz  = 

J     O      m      A*2  J  n    >  n       h    n 


2K  Vhkn  6^  -  {AkVhk  Zn  Zm  dz 


(28) 


In  order  to  recast  the  first  term  of  Eq.  (28)  in  a  form  useful 
for  this  derivation,  we  first  use  integration  by  parts  twice 
with  respect  to  z.  The  resulting  expression,  after  some 
lengthy  manipulations,  is 


CO 

/ 


lz  d2VA 
P    m    dz2 


*°  -  u  ^ d2z- 


dz: 


dz 


Pi    lm      dz         0 


|ff(r,6) 


j^z   dvAz2n  _   i  az 

P2    2"      5^        H(r,0)        Pi    3z 


lm  y   z 


\H(r,Q) 
0 


1    dz2/nv7 


i/(r,6) 


(29) 
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Again,  the  subscripts  1  and  2  denote  the  water  column  and  the 
sediment,  respectively. 

Application  of  the  boundary  conditions  Eqs .  (6),  (7),  (8), 
(14),  (15),  (17),  (18)  and  (19)  to  Eq.  (29),  yields 
subsequently 


/- 

J     O 


1  „   d\2Zn 
ozz 


dz  = 


/>a— *  <*■ ♦ 


0 


dz: 


J_z    {dv^ln    dvhz2n] 


'irn 


dz 


dz 


z=H(r,Q) 


dz 


IB  \  J^XJ  Z     -  — V  Z     ' 

dz    L  9l   h  ln    P2  h  2n 


z=H(r,Q) 


(30) 


Replacing  the  first  term  of  Eq.  (28)  by  Eq.  (30)  and  then 
making  use  of  the  depth  equation  (Eq.  (3)),  we  obtain  the 
following  alternative  expression  for  the  vector  mode  coupling 
coefficient,  for  iw*-n    : 


15 


A  A 


A7      ,M5      d\Zzn, 
9lXmK      dz  dz     ' 


knHr,Q)-km2(r,e) 


[2Iik  v^ z» z-  dz  + 


z=Hlr,Q) 


dZi"  flu  Z     -  — V  Z     ) 

az  (  Pl  *  ln   p2  ^  2n) 


z=H(r,Q) 


(31) 


or  equivalently ,  in  light  of  the  boundary  conditions  Eq.  (14) 
and  Eq.  (15) , 


A«(*,e)  -  P^r  +  Ymne  = 


*aa(r,-6)-V(^e) 


[2Iik  Whk  z« z"  dz  + 


-L(l-fi)Z 

Pi  P2 


=//(r,0)         ( 


i     i  ,  a^.  VkZ 


Pi      ?2       dz 


h^m 


--H(x.Q)    -I 


(32) 


The  above  expression  only  involves  Zn  and  not  their  horizontal 
derivatives  in  the  integrands.  Therefore,  the  corresponding 
numerical  evaluations  are  more  efficient. 

The  last  two  terms  of  Eq.   (32)   express  the  direct 
contribution  of  bathymetry  change  and  sediment  properties  in 

— * 

A^  .  They  were  excluded  in  the  previous  model  but  are  included 
in  the  latest  version. 
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Jb.  Case  of  m=n 

In  order  to  derive  an  expression  for  the  vector 
mode  coupling  coefficient  for  m=n,  we  differentiate  the 
orthonormal  condition  Eq.  (4)  using  Leibniz  rule.  The  result 
is 


A     A 
AM(rf0)  =  pmr  +  Ynn8  = 


(33) 


-  VJf(r,6)  (  — -  — )  z 


±-±)  z   2 

Px  P2   ln  'z=H(r,d) 


Note  that,  this  coupling  coefficient  is  zero  for  a  flat 
horizontal  bottom.  The  latest  version  of  the  model  has 
included  this  new  term. 

2.   SCALAR  MODE  COUPLING  COEFFICIENT,  Bmn 

Taking  the  horizontal  gradient  of  both  sides  of  Eq. 
(21),  i.e.,  definition  of  the  vector  mode  coupling 
coefficient,  and  applying  the  Leibniz  rule  for 
differentiation  of  a  definite  integral,  we  get  after  some 
manipulations,  the  following  expression  for  the  scalar  mode 
coupling  coefficient: 
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B*n<r,e)  =  jvL(r,6)    f±  VhZm .  V hZ  n  dz   - 

Z  o  P  (34) 

Pi  P2   lin  h  ln     ^  z=H(r,d) 


There  is  a  unique  property  associated  with  a  complete  set  of 
orthonormal  functions,  called  the  "closure  relationship."   For 
the   normal  modes,  this  relationship  can  be  expressed  as 

£  ±Zn(z;r,e)    Zn(z/;rf6)  =b{z-z')  <35> 

n     p  n  n 


Taking  the  horizontal  gradient  of  both  sides  of  the  closure 
relationship,  multiplying  by  Zm(z;r,0)  and  then  integrating 
over  the  entire  depth,  we  get,  after  some  rearranging  of 
terms, 

VhZm(z;r,Q)    =   -  E  E^  (r,  6)  Zn  (z;  r,  6)  (36, 


where 
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Em{r,0)    =  ^A^U.e)    + 


VhH(r,d)  (±-±)Zln(H;r,Q)Zlm(H;r,d) 
Pi  P2 


(37) 


Substituting  now  Eq.  (36)  in  Eq.  (34),  we  finally  obtain  the 
following  alternative  expression  for  B^: 


B«  =  ^V^^Cr.e)  -EEnJ(r,6).i?mJ(r,e) 

(38) 


(  —  -—  )Zljn(z;r#e)V/IZln(z;r#e)  .VAff(r,8; 
Pi   P2 


z=H(r,Q) 


Eq.  (38)  is  valid  for  both  the  rrp^n  and  m=n  cases.  The  last 
term  of  Eq.  (38),  is  new  in  the  model.  The  magnitude  of  this 
mode  coupling  coefficient  is  generally  much  less  than  the 
magnitude  of  the  vector  coefficient. 

C.   NUMERICAL  IMPLEMENTATION 

The  major  part  of  the  model  upgrade  was  the  replacement  of 
the  old  routines  with  new  ones  for  the  computations  of  the 
rederived  mode  coupling  coefficients  according  to  Eqs .  (32), 
(33)  and  (38) . 

These  new  routines  are  contained  in  a  program  called 
"sedbot"  and  are  listed  in  Appendix  C.  Normal  modes  and  the 
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horizontal  gradient  vector  of  the  wavenumber  as  function  of 
position  are  required  as  input  to  "sedbot. " 

The  normal  mode  field  is  created  by  a  program  called 
"modes",  whereas  the  horizontal  gradients  of  the  wavenumber 
are  calculated  in  the  program  "kder. "  These  two  programs  are 
listed  in  Appendices  A  and  B,  respectively. 

D.   EXAMPLE  RUNS 

For  both  example  runs,  the  medium  is  taken  to  have  two 
isospeed,  isodensity  layers  separated  by  a  constant -slope 
interface.  Sound  speed  is  taken  to  be  1500  m/sec  in  the  water 
column  and  3000  m/sec  in  the  sediment.  Density  is  taken  to  be 
1000  Kg/m3  in  the  water  column  and  2500  Kg/m3  in  the  sediment. 
The  source  is  taken  to  be  harmonic  in  time,  with  a  frequency 
of  100  Hz,  and  is  located  at  a  depth  of  50  m.  The  coupled  mode 
solution  along  two  radii,  90°  and  45°,  are  displayed  and 
discussed. 

1.   BOTTOM  SLOPE  =  .001  RADIANS 

The  bottom  slope  in  this  first  case  is  taken  to  be 
.001  radians.  The  water  depth  is  100  m  at  the  source  location 
and  70  m  after  30  km  in  the  y  direction  (see  Fig.  2)  . 

At  the  source  location  there  are  twelve  trapped  modes  in 
the  water  layer.  Only  eight  trapped  modes  exist  at  the 
location  30  km  upslope. 
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p-1000Kg/m   '  p  =  2500  Kg/m 

w  s 

c  =  1 500  m/sec  c  =  3000  m/sec 

w  s 


a) 


100  m 


pressure  release  surface        30  Km 


WATER  LAYER 


/|\ 

70  m 

V 


SEDIMENT  LAYER 


X  / 

15  Km 


b) 


source 


-15  Km 


30  Km    y 


Figure  2.    Geometry  of  the  first  example  case  with  a 
constant  slope  of  .001  radians  along  y-axis  (a),  and  a  plane 
view  showing  the  6    =90°  and  6    =  45°  propagation  paths  (b) 
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Figure  3.    Envelope  amplitudes  of  the  first  eight  trapped 
modes  in  the  3D  coupled  mode  solution  along  the  path  0    =  90° 
for  a  bottom  slope  of  .001  radians 
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Figure  4.  Envelope  phases  of  the  first  eight  trapped  modes 
in  the  3D  coupled  mode  solution  along  the  path  $  =90°  for  a 
bottom  slope  of  .001  radians 
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Figure  5.    Envelope  amplitudes  of  the  first  eight  trapped 
modes  in  the  3D  coupled  mode  solution  along  the  path  0    =45 
for  a  bottom  slope  of  .001  radians 
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Figure  6.  Envelope  phases  of  the  first  eight  trapped  modes 
in  the  3D  coupled  mode  solution  along  the  path  6  =  45°  for  a 
bottom  slope  of  .001  radians 
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Fig.  3  and  Fig.  4  show  the  amplitude  and  phase  of  the 
modulation  envelope,  Un,  for  the  first  eight  modes  travelling 
in  the  upslope  y  direction,  i.e.,  along  the  path  6  =  90°  (see 
Fig.  2 )  .  An  upslope  enhancement  is  noticed,  especially  for  the 
higher  order  modes,  as  they  propagate  into  shallower  water. 
The  phase  of  the  envelope,  which  is  the  phase  deviation  from 
the  2D  adiabatic  approximation,  is  very  small  (about  11° 
maximum) .  The  amplitude  fluctuations  are  between  15%  and  30% 
for  all  the  modes.  In  light  of  the  small  amplitude  and  phase 
fluctuations,  the  adiabatic  approximation  can  be  considered 
reasonable  along  this  propagation  path. 

Fig.  5  and  Fig.  6  show  the  amplitude  and  phase  of  the 
modulation  envelope,  Un,  for  the  first  eight  modes,  along  the 
propagation  path  6  =  45°  (see  Fig.  2)  .  Here,  the  upslope 
enhancement  is  significantly  less  and  the  fluctuations  of  the 
amplitude  of  the  higher  order  modes  at  greater  range  are 
slightly  larger  than  along  the  previous  path.  We  speculate 
that  this  slight  increase  in  the  fluctuations  is  due  to  that 
more  interacting  modes  remain  trapped  in  the  water  column  at 
longer  ranges  along  this  path.  The  higher  order  modes  have 
large  phase  deviations  from  the  2D  adiabatic  phases.  These 
large  phase  changes  correspond  to  significant  horizontal 
refraction  of  the  wave  fronts  due  to  the  existence  of  a 
transverse  gradient  in   the  bottom  bathymetry. 
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Figure  7.    Geometry  of  the  second  example  case  with  a 
constant  slope  of  .002  radians  along  y-axis  (a),  and  a  plane 
view  showing  the  $    =90°  and  9    =  45°  propagation  paths  (b) 
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Figure  8.    Envelope  amplitudes  of  the  first  eight  trapped 
modes  in  the  3D  coupled  mode  solution  along  the  path  0  =  9  0 
for  a  bottom  slope  of  .002  radians 
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Figure  9.  Envelope  phases  of  the  first  eight  trapped  modes 
in  the  3D  coupled  mode  solution  along  the  path  $  =90°  for  a 
bottom  slope  of  .002  radians 
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Figure  10.    Envelope  amplitudes  of  the  first  eight  trapped 
modes  in  the  3D  coupled  mode  solution  along  the  path  6    =45° 
for  a  bottom  slope  of  .002  radians 
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Figure  11.    Envelope  phases  of  the  first  eight  trapped 
modes  in  the  3D  coupled  mode  solution  along  the  path  6    =45 
for  a  bottom  slope  of  .002  radians 
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2.   BOTTOM  SLOPE  =  .002  RADIANS 

For  this  case,  the  same  isospeed,  isodensity,  wedge 
shape  waveguide  is  used,  except  the  bottom  slope  is  now 
doubled  (.002  radians).  The  bottom  depth  at  the  source 
location  is  now  100  m  and  shoals  to  40  m  after  30  km  away  from 
the  source  in  the  y  direction  (see  Fig.  7)  .  At  the  source 
position  there  are  twelve  trapped  modes,  but  in  30  km  upslope, 
there  are  only  five  trapped  modes  in  the  water  column. 

Fig.  8  and  Fig.  9  show  the  amplitudes  and  phases  of  the 
modulation  envelope,  Un,  for  the  first  eight  trapped  modes 
travelling  upslope  in  the  y-axis  direction,  i.e.  along  the 
path  9  =  90°  (see  Fig.  9)  .  Upslope  enhancement  is  much 
stronger  than  the  previous  case,  especially  for  the  higher 
modes.  The  fluctuations  in  amplitude  is  about  50%  in  some 
modes  and  in  phase  more  than  20°.  Thus,  the  adiabatic 
approximation  would  induce  considerably  larger  errors  than  the 
case  of  a  .001  bottom  slope. 

Fig.  10  and  Fig.  11  show  the  amplitudes  and  phases  of  the 
modulation  envelope,  Un,  for  the  same  eight  modes  along  the 
propagation  path,  i.e.,  6  =  45°  (see  Fig.  7).  The  horizontal 
refraction  phenomenon  is  much  stronger  here  than  for  the  case 
of  a  .001  slope.  Along  this  path,  the  adiabatic  approximation 
would  also  induce  large  errors.  Typical  percentages  of 
amplitude  fluctuations  are  about  50%  for  the  second  mode  and 
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30%  for  the  third  and  fourth  modes.  The  phase  deviation, 
especially  for  the  higher  order  modes,  is  also  large. 
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IV.    CONCLUSIONS  AND  RECOMMENDATIONS 

A  3D  coupled  normal  mode  model  for  sound  propagation  in 
shallow  water  with  irregular  bottom  bathymetry,  is  developed 
in  this  thesis.  This  model  can  be  used  to  examine  underwater 
sound  propagation  involving  significant  bottom  interaction.  In 
this  model,  sound  speed  is  allow  to  vary  in  three  dimensions 
and  water  depth  and  sediment  properties  in  horizontal 
location . 

It  is  shown  here  that,  for  a  frequency  of  100  Hz,  the 
adiabatic  approximation  is  valid  only  for  very  mild  bottom 
slopes.  Typical  errors  for  a  slope  of  .001  radians  are  15%  in 
mode  amplitude  and  10°  in  its  phase.  For  a  slope  of  .002 
radians,  the  errors  are  significantly  larger. 

The  model  presented  in  this  thesis  is  capable  of 
simulating  the  interactions  of  the  normal  modes  as  they 
propagate  in  complex  environments.  Propagation  phenomena  like 
mode-mode  interaction,  horizontal  refraction  and  slope 
enhancement  can  be  examined  using  this  improved  model. 

In  the  development  of  the  present  model  an  approximation 
(Eq.  12)  in  the  bottom  boundary  conditions  is  used.  The 
validity  of  this  approximation  requires  that  the  slope  must  be 
much  smaller  than  unity  (Eq.  13)  .  In  order  to  be  able  to 
handle  very  steep  bottom  slopes,  i.e.,  order  one  slopes,  one 
needs  to  use  the  exact  form  of  the  bottom  boundary  conditions 
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(Eq.  10)  .  This  could  make  the  formulation  more  complicated  but 
it  should  be  tractable. 

Another  future  improvement  to  the  model  will  be  to  include 
sound  energy  absorption  (attenuation) .  One  way  to  do  this  is 
by  introducing  imaginary  parts  in  the  eigenvalues 
(wavenumbers)  .  Lastly,  a  test  for  the  accuracy  of  the  improved 
model  is  needed.  This  can  be  achieved  by  comparing  the  results 
generated  by  this  model  with  some  exact  analytic  solutions. 
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APPENDIX  A.   FORTRAN  ROUTINES  FOR  COMPUTING  NORMAL  MODES 
FIELD 

The  following  program  creates  the  normal  mode  field, 
"amode.dat",  for  a  given  geographical  area  of  the  ocean.  Sound 
speed,  density,  and  bottom  depth,  are  defined  for  every  grid 
point.  Given  sound  speed  field  and  density,  the  normal  modes 
are  calculated  using  a  standard  mode  solver  routine. 

************************************************************ 

*  * 

*  This  program  computes  the   normal  modes  in  a   given     * 

*  area.  * 

*  • 
******************  ****************************************** 

c 

C  INPUT/ARGUMENTS 

c  f  frequency,  Hz 

c  xmax  maximum  position  in  x  direction,  meters 

c  ymax  maximum  position  in  y  direction,  meters 

c  h(nx,ny)  depth  at  x,y  position,  meters 

c  nx  number  of  stations  in  x  direction 

c  ny  number  of  stations  in  y  direction 

c  nz  number  of  stations  in  z  direction 

c  xo  initial  x  position,  meters 

c  yo  initial  y  position,  meters 

c  dz  step  size  in  z  direction,  meters 

c  cf (nx,ny,nz)  sound  speed  field  in  every  x,y  position, 

c  m/s 

c  df (nx,ny,nz)  density  field  in  every  x,y  position, 

c  kgr/m^3 

c  isw  switch  index  :  1  write  /  0  do  not  write 

c  mm  maximum  number  of  allowable  trapped 

c  modes 

c 

c  OUTPUT/ARGUMENTS 

c  for  each  horizontal  station  x,y: 

c  ksq_r(nz)  squared  eigenvalues  for  each  trapping 

c  mode  (real) 

c  efun_r (nz , nz)  eigenfunctions  (real) 

c  h(nx,ny)  water  depth  (meters) 

c  c(nz)  sound  speed  profile  in  a  specific  grid 
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c  position 

c  d(nz)            density  profile  in  a  specific  grid 

c  position 

c  w               source  angular  frequency  (rad/sec) 

c  dx,dy,dz        step  size  in  x,y  and  z  directions  (meters) 

c  nx,ny,nz        number  of  stations  in  x,y  and  z  directions 

c  (metrers) 


c 


c 

c 


program  modes 

parameter (xmax=3  0000 . dO , ymax=3  0000 . dO , nx=ll, ny=ll , 

*  nz=100,mm=20) 
implicit  real*8  (a-h,o-z) 

real* 8  cf (nx, ny, nz) , c (nz) , df (nx, ny, nz) , d (nz) 
real*8  h (nx, ny) , ksq_r (nz) , ksq_i (nz) 
real* 8  efun_r (nz, nz) , efun_i (nz, nz) , efun (nz,mm) 
real* 8  ks (mm) , x (nx) , y (ny) , z (nz) , ksed, kwat 
logical  ex 

data  isw  /l/ 

inquire ( f ile= ' amode . dat ' , exist=ex) 
if  (ex)  then 

open (unit =13 , f ile= ' amode .dat ' , status=' old' ) 

close (13 , status=' delete' ) 
endif 
open (unit=4 , f ile= ' /home/noise/sagos/modes/amode . dat ' , 

*  form=' unformatted' , status^ ' new' ) 
inquire (f ile='mode. sys ' , exist=ex) 

if  (ex)  then 

open (unit =13 , f ile='mode . sys' , status= 'old' ) 

close ( 13 , status = 'delete' ) 
endif 
open (unit=6 , f ile= ' /home/noise/sagos/modes/mode . sys ' , 

*  form= ' formatted' , status= ' new' ) 

write (6 ,*)' output  field' 

dx=xmax/df loat (nx- 1) 

dy=ymax/df loat (ny- 1) 

dz=2.d0 

pi=4 . dO*datan ( 1 . dO ) 

f=224.d0 

w=2 .dO*pi*f 

call  data  (cf ,  df ,  h,  nx,  ny,  nz,  dx,  dy,  dz,  xo,  yo,  x,  y,  z) 

nm=nz-2 

write(4)  w,dx,dy,dz 

write(4)  nx,ny,nz 

if  (isw.eq.l)  then 
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write (6 , *) ' dz=  ',dz,  '  meters' 
write (6 , *) ' nx=  ',nx,'  ny=  ',ny,'  nz  =  ' , nz 
write (6 ,*)' interface  depth' 
write (6, *) ' h=  ',h(l,l)  ,'  meters' 
write (6 ,*)' sound  speed  profile,  m/s' 
write (6 , *) (cf(l,l,iz) ,iz=l,nz,2) 
write (6 ,*)' density  profile,  kgr/m^3' 
write (6 , *) (df (1, 1, iz) ,iz=l,nz,2) 
endif 
c 

mcntr=0 
c 

do  ix=l,nx 

do  iy=l,ny 
icounter=0 
ichk=0 

if  ( ix. eq. 1 .and. iy . eq. 1)  ichk=4 
write (4)  h(ix,iy) 
do  iz=l,nz 

c (iz) =cf (ix, iy, iz) 
d (iz) =df (ix, iy, iz) 
enddo 

write(4)  c 
write(4)  d 

write (6 , *) ' ix=' , ix, '  iy=',iy,'  ichk=',ichk 
call  mode ( f , nz , dz , c , d, nm, ksq_r , ksq_i , 

*  efun_r, efun_i, ichk) 
c 

c    choose  only  the  trapped  modes 

c    cs  :  sound  speed  in  the  sediment  (constant) 

c    cw  :  sound  spedd  in  the  water  column, next  to  the 

c        interface 

c 

cs=c (int (h (ix, iy) /dz) +2) 
cw=c (int (h (ix, iy) /dz) - 1) 
c 

c    set  zeros  in  the  eigenvalues -eigenf unctions  arrays 
c 

do  i=l,mm 

ks  (i)=0.dO 
do  j=l,nz 

efun( j , i) =0.d0 
enddo 
enddo 
c 

do  i=l,nm 

ksed= (w/cs) **2 
kwat= (w/cw) **2 
if  (ksed. It .ksq_r (i) .and.ksq_r (i) .lt.kwat) 

*  then 

icounter=i count er+l 
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ks ( icounter) =ksq_r ( i) 
do  iz=l,nz 

ef un ( iz , icounter) =efun_r ( iz , i) 
enddo 
endif 
enddo 
c 

if  (isw.eq.l)  then 

write(6, *) ' h=  ',h(ix,iy) 

write (6 ,*)' limits  for  trapped  modes  :  ' 
write ( 6 , * ) kwat , ksed 
write (6 , *) ' icounter=  ', icounter 
write (6, *) 'kn^2' 
write(6,*) (ks(i) ,i=l,mm) 
endif 
c 

if  (ix. eq. 3 .and. iy . eq. 3 )  then 
write (6, *) ' efun(iz, 18) ' 
write(6,*) (efun(iz,18) ,iz=l,nz) 
endif 
c 

write (4)  icounter 

if  (icounter . gt .mm. and. icounter . gt .mcntr) 
*  mcntr=icounter 

write (4)  ks 
write (4)  efun 
enddo 
enddo 
c 

if  (mcntr. ne.O)  write (6 , 1001) mcntr 
1001  format (i3,'  trapped  modes,  exceeds  limit,  increase  mm 
c  and  rerun' ) 

close (4) 
close (6) 
end 


************************************************************ 

*  * 

*  The  following  program  provides  an  example  data  input.    * 

*  * 

************************************************************ 


c 

c 

IN 

c 

nx 

c 

ny 

c 

nz 

c 

dx 

c 

dy 

c 

dz 

INPUT/ARGUMENTS 


number  of  stations  in  x  direction 
number  of  stations  in  y  direction 
number  of  stations  in  z  direction 
step  size  in  x  direction,  meters 
step  size  in  y  direction,  meters 
step  size  in  z  direction,  meters 
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c    xo  initial  x  position,  meters 

c    yo  initial  y  position,  meters 

c 

c    OUTPUT/ARGUMENTS 

c    cf (nx,ny,nz)     sound  speed  field  in  every  x,y  position, 

c  m/s 

c    df (nx,ny,nz)     density  field  in  every  x,y  position, 

c  kgr/mA3 

c   h(nx,ny)        interface  depth,  meters 

c 

subroutine  data  (cf  ,  df  ,  h,  nx,  ny,  nz,  dx,  dy,  dz,  xo,  yo,  x,  y ,  z) 
c 

implicit  real*8  (a-h,o-z) 

real* 8  cf (nx, ny, nz) , df (nx, ny, nz) , h (nx, ny) 

real*8  x (nx) , y (ny) , z (nz) 
c 

do  ix=l,nx 
do  iy=l,ny 

x (ix)  =xo+dx*df loat (ix- 1) 
y (iy) =yo+dy*df loat (iy- 1) 

c 

c    bathymetry  field 

c 

h(ix, iy) =- .0  005d0*x(ix) +10  0.d0 
c 

c    sound  speed  and  density  fields 
c 

do  iz=l,nz 

z (iz) =dz*dfloat (iz-1) 
if  (z  (iz)  . le.h (ix, iy) )  then 
cf (ix, iy, iz) = . 0005d0*x ( ix) - 
*  .ldO*z (iz) +1490. dO 

df (ix, iy, iz) =1000.d0 
else 

cf(ix,iy,iz)=1800.d0 
df (ix, iy, iz) =2  000.d0 
endif 
enddo 
enddo 
enddo 
c 

return 
end 
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APPENDIX  B.   FORTRAN  ROUTINE  FOR  COMPUTING  WAVENUMBER 
DERIVATIVES 

This  program  inputs  from  "amode.dat"  as  created  by  the 

previous  program  "modes" .   Computes  the  horizontal  derivatives 

of  the  total  wavenumber  k,  at  every  position  of  the  acoustic 

field.  The  derivative  calculation  requires  definition  of  a 

computational  domain.  Output  is  to  the  file  "kder . dat " . 

************************************************************ 

*  * 

*  This  program  assigns  the  source  position  relative  to      * 

*  the  input  field  via  ixorig  and  iyorig.  * 

*  This  program  also  specifies  the  radial  increment  for  the  * 

*  spline  definition,  the  number  of  intervals,  and  the      * 

*  angular  increment   between  integration  paths  (dr,nr,da)   * 

*  * 

*  Procedure:  xy- spline  at  each  depth  * 

*  evaluate  dk/dx,dk/dy  * 

*  transform  into  dk/dr,dk/da  * 

*  * 
******************************************************** 

c 

program  kder 
c 

parameter (nx=ll, ny=ll , nz=100 , nm=20 , ndum=nz*nm, 
*         ixorig=l , iyorig=3 , nwk=2*ny*nx+2*max0 (nx, ny) ) 

implicit  real*8  (a-h,o-z) 

real *  8  kd ( nx , ny , nz , 2 )  , k ( nx , ny ) , kdxy ( 6 ) , kc ( 2 , nx , 2 , ny ) 

real* 8  x(nx) , y (ny) , ang (nx, ny) , c (nz) , ct (nx, ny, nz) 

real* 8  wk(nwk) ,efun(ndum) ,hork(nz) 

character*20  filename 

logical  ex 
c 
c 

c  open  statements 

c 

c 

inquire (f ile=' kder .dat ' , exist =ex) 
if  (ex)  then 

open (unit =13, file=' kder .dat ' , status=' old' ) 

close (13 , status=' delete' ) 
endif 
inquire (f ile=' kder .sys' , exist=ex) 
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if  (ex)  then 

open (unit =13 , FILE=' kder . sys ' , status=' old' ) 

close (13 , status=' delete' ) 
endif 

inquire ( f ile= ' ef un_orig . dat ' ; exist=ex) 
if  (ex)  then 

open (unit=13 , f ile=' efun_orig .dat' , status=' old' ) 

close (13 , status = ' delete' ) 
endif 

open (unit=4 , f ile= ' amode . dat ' , f orm= ' unformatted' , 

*  status='old' ,err=2001) 

open (unit=14 , f ile= ' kder . dat ' , f  orm= ' unformatted' , 

*  status='new' ,err=2002) 

open (unit=24 , f ile= ' ef un_orig . dat ' , f orm= 

*  'unformatted' , status= ' new' ,err=2003) 

open (unit =6 , f ile=' kder. sys' , f orm=' formatted' , 

*  status=' new' , err=2004) 
c 

c  input  w  (rad/sec) ,  dx,dy  (meters),  dz  (meters) 

read (4)  w,dx,dy,dz 
c 

c  input:  number  of  x  indices,  no.  y  indices 
c        number  of  modes,  TOTAL  vertical  increments 

read (4)  nxt,nyt,nzt 
c 

write (6, 1009)  w 

write ( 6 , * ) '  dx (m) , dy (m) , dz (m) ' 

write(6,*)  dx,dy,dz 

if  (m.ne.nm)  write  (6  ,  *)  'm='  ,m,  '   nm=',nm 

if  (nx.ne.nxt)  write (6 , *) ' nx=' , nx, '   nxt=',nxt 

if  (ny.ne.nyt)  write (6 , *) ' ny=' , ny, '   nyt=',nyt 

if  (nzpl .ne.nzt)  write (6 , *) ' nzpl= ' , nzpl, '   nzt=',nzt 
c 

pi  =  dacos(-l.dO) 
c 
c  distribution  parameters  for  spline 

dda  =  .8d0 

da  -  dda*pi/180. 

dr  =  3000. dO 
c  number  of  points  in  spline 

nr  =  (nx-ixorig) *dx/dr  +1 

write (6,*)  nr, '  spline  locations  with  interval=' , dr 
c  number  of  radial  paths 

tang  =  datan (  .5*dy*(ny-l)  /  (dx*(nx-l))  ) 
c      na  =  2*idint (tang/da)  +  1 

na=3 

write (6,*)  na, '  integration  paths  for  da=',dda,'  deg' 
c 

c 

c  horizontal  field  grid  in  meters 
do  11  ix=l,nx 
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11  x(ix)  =  df loat (ix- ixorig) *dx 
do  12  iy=l,ny 

12  y(iy)  =  df loat (iy- iyorig) *dy 
write(6,*)'x  range:  ( ' ,x (1) , ' , ' ,x (nx) , ' ) ' 
write(6,*)'y  range:  ( ' , y (1) , ' , ' , y (ny) , ' ) ' 

c 

c  read  in  c- field  and 

c  calculate  angle  (ccw  from  x-axis) 

do  14  ix=l,nx 

do  14  iy=l,ny 
c  full  sound  speed  prof ile (nzpl)  (0-5000m) 

read (4)  h 

read (4)  c 

read (4)  dens 

read (4)  icounter 

read (4)  hork 

read (4)  efun 
c  create  file  to  obtain  initial  conditions 

if  (iy. eq. iyorig. and. ix. eq. ixorig)  write (24) efun 

do  iz=l,nz 

ct(ix,iy,iz)  =  c(iz) 

enddo 

if  (ix. eq. ixorig)  then 

if  (iy .ge . iyorig)  ang(ix,iy)  =  pi/2. 

if  (iy. It . iyorig)  ang(ix,iy)  =  -pi/2. 

else 

ang(ix,iy)  =  datan ( (y (iy) ) / (x (ix) ) ) 

endif 
14  continue 
c 

c  calculate  derivatives  from  first  station  below  surface  to 
c  bottom 
c 

do  100  iz=l,nz 

do  110  ix=l,nx 
do  110  iy=l,ny 
c 

c  wavenumber  is  in  rad/m 
c 

110     k(ix,iy)  =  w/ct (ix, iy, iz) 
c 

c  fit  bi-cubic  spline  to  iz-th  level  waveno. 
c 

ic  =  nx 

call  ibcccu (k, x, nx, y, ny, kc, ic, wk, ier) 

if  (ier.ne.0)  write (6 , 1001)  ier 
c 

c  use  spline  to  evaluate  cartesian  derivatives  at  each  grid 
c  point 

c  transform  derivatives  into  cylindrical  coordinates 
c 


43 


do  12  0  ix=l,nx 
do  121  iy=l,ny 

call  dbcevl (x, nx, y, ny, kc , ic, x (ix) , y (iy) , kdxy , ier; 
if  ( ier . ne . 0 ) 

*  write (6, 1007)  ix, iy , X ( ix) , Y ( iy) , ier 
call  cyl (ang ( ix, iy) , kdxy (2 ) , kdxy (3 ) , 

*  kd (ix, iy, iz, 1) , kd (ix, iy, iz, 2 ) ) 
121       continue 

120    continue 

if  (iz.eq.l)  then 

write (6,*)'  For  first  subsurface  layer' 

write (6, 1005) 

write (6,1006)  ( (ct(ix,iy,l) , ix=l, nx) , iy=ny, 1,-1) 

write (6, 1011) 

write (6,1016)  ( (k (ix, iy) , ix=l, nx) , iy=ny, 1,-1) 

write (6, 1003) 

write (6, 1002)  ( (kd (ix, iy , 1, 1) , ix=l, nx) , iy=ny, 1,-1) 

write (6, 1004) 

write (6, 1002)  ( (kd (ix, iy, 1, 2) , ix=l,nx) , iy=ny, 1,-1) 

write (6, *) 
endif 


100  continue 


write 
write 
write 
write 

write 
write 
write 
write 
write 
write 
write 
write 
write 
write 
write 


14)  ixorig, iyorig 

14)  da,na,dr,nr 

6,*)  ' da, na, dr, nr' , da, na, dr, nr 

14)  kd 

6,*)'  For  bottom  level' 

6,1005) 

6,1006)  ( (ct (ix, iy, nz) , ix=l, nx) , iy=ny, 1,-1) 

6, 1011) 

6, 1016)  ( (k(ix, iy) , ix=l,nx) , iy=ny, 1,-1) 

6. 1003) 

6, 1002)  ( (kd(ix, iy,nz, 1) , ix=l,nx) , iy=ny, 1,-1) 

6. 1004) 

6, 1002)  ( (kd(ix, iy,nz, 2) , ix=l,nx) , iy=ny, 1,-1) 

6,*) 

6, 1008)   (kd(4,3,iz,l) ,iz=l,nz,2) ,kd(4,3,nz,l) 


goto  2020 

1001  f ormat ( '  ier:',i3,'  for  ibcccu,  xy-spline') 

1002  format (5 (lx,el2.5) ) 

1003  formate  dk/dr   (  (rad/m)   /  m  )') 

1004  formate  dk/rda   (  (rad/m)   /   m  )') 

1005  format ( '  c   (m/s) ' ) 

1006  format (5 (2x, f 8.3)  ) 
1016  format (5 (2x, f8.5)  ) 

1007  formate  ix,  iy,X,  Y,  ier  for  dbcevl:  '  ,  213  ,  2F6  . 1,  i3) 

1008  formate  at  ix,iy=4,3    dk/dr  (z)  ' /ll  (5  (lx,  dll .  5) /) 
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1009  formatC  frequency,  rad/sec  :',dl4.7) 

1010  formatC  k'  /7  (4  (3x,  dll .  4)  /3x,  3  (3x,  dll .  4)  / 

1011  formatC  k   (rad/m)  '  ) 

c 

c 

c  close  statements 

c 

c 

2001  f ilename=' amode . dat ' 
goto  2010 

2002  f ilename=' kder.dat ' 
goto  2010 

2003  filename=' efun_orig.dat ' 
goto  2010 

2004  f ilename=' kder . sys' 

2010  write (*, 2011)  filename 

2011  formatC  ERROR  OPENING  FILE  :  ',A) 
2020   close(4) 

close (14) 
close (6) 
end 


subroutine  cyl (ang,x,y, r, a) 
polar  transformation  subroutine 
implicit  real*8  (a-h,o-z) 
r  =  x*dcos(ang)  +  y*dsin(ang) 
a  =  x*dsin(ang)  -  y*dcos (ang) 
return 
end 
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APPENDIX  C. 
COEFFICIENTS 


FORTRAN  ROUTINES  FOR  COMPUTING  THE  COUPLING 


************************************************************ 


* 
f\* 


This  program  manages  the  subroutines  "subl.f",  "sub2 

"sub3.f",  "partial. f"  which  compute  the  two  mode  * 

coupling  coefficients.  The  input  is  from  "amode.dat"  * 

and  "kder .dat " ,  specifically  the  modes,  horizontal,  * 

wavenumber,  horizontal  derivatives  of  total  wavenumber,  * 
bathymetry  and  density.  The  output  file  is  "mcoupl.dat".* 

*  * 

************************************************************ 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


rhol  :  water  column  density- constant  in  depth  (kg/m^3) 

rho2  :  sediment  density- constant  in  depth  (kg/m^3) 

dz  :  vertical  step  size  (m) 

run  :  maximum  number  of  trapped  modes, in  the  water  column 

nx,ny  :  number  of  stations  in  x  and  y  directions 

h   :  bottom  bathymetry 

zb  :  acoustic  pressure  eigenf unctions, at  the  interface 

depth 


zbml 

zbm  : 

zbmml 

zbn  : 

zbnml 

ar 

aa 

cr 

ca 


acous 

zb  for 
:  zbml 
zb  for 
:  zbml 
range  c 
angle  c 
correct 

correct 


b   :  second 
k   :  square 


tic  pressure  eigenf uctions, one  step  size 

above  the  interface  depth 
the  mth  mode 
for  the  mth  mode 
the  nth  mode 
for  the  nth  mode 
omponent  of  the  first  coupling  coeff. 
omponent  of  the  first  coupling  coeff. 
ion  at  the  range  component  of  the  first 

coupling  coeff. 
ion  at  the  angle  component  of  the  first 

coupling  coeff. 
coupling  coeff. 
of  horizontal  wavenumbers (eigenvalues)  of 

the  modes 


program  sedbot 

parameter (nx=5 , ny=5 , nz=100 , nw=2*nx*ny+2*max (nx, ny) , 
*  nm=2  0) 

implicit  real*8  (a-h,o-z) 

real* 8  h  (nx,  ny)  ,  zbm(nx,  ny)  ,b  (run,  nm,  nx,  ny)  ,  km(nx,  ny) 
real*8  zbmml  (nx,  ny)  ,  cr  (nm,  nm,  nx,  ny)  ,  ca  (nm,  nm,  nx,  ny) 
real*8  ar (nm, nm, nx, ny) , aa (nm, nm, nx, ny) ,kn(nx,ny) 
real* 8  zbn (nx, ny) , zbnml (nx, ny) , c (2 , nx, 2 , ny) , wk (nw) 
real* 8  crl (nx, ny) , cal (nx, ny) , x (nx) , y (ny) 
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real *  8  zbnpr ( nx , ny )  , zbnpa ( nx , ny )  , zb ( nm , nx , ny ) 

real *  8  zbnpx ( nx , ny )  , zbnpy ( nx , ny )  , hpr ( nx , ny )  , hpa ( nx , ny ) 

real*8  hpx (nx, ny) , hpy (nx, ny) ,bl (nx, ny) , kk (nm) 

real* 8  rhol (nx, ny) , rho2 (nx, ny) , k (nx, ny, nz) 

real* 8  cs (nz) , d (nz) , efun (nz , nm) , kd (nx, ny, nz , 2 ) 

integer  icounter (nx, ny) 

logical  ex 

inquire (f ile='mcoupl . dat ' , exist=ex) 
if  (ex)  then 

open (unit =13 , f ile= 'mcoupl .dat ' , status= 'old' ) 

close (13 , status = ' delete' ) 
endif 
open (unit=4 , f ile= ' /home/noise/sagos/modes/amode . dat ' , 

*  form= ' unformatted' , status=' old' ) 
inquire ( f ile= ' coupl . sys ' , exist=ex) 
if  (ex)  then 

open (unit=13 , f ile= ' coupl . sys ' , status= ' old' ) 

close (13 , status=' delete' ) 
endif 
open (unit=6 , f ile= ' /home/noise/sagos/modes/coupl . sys ' , 

*  form=' formatted' , status=' new' ) 
open (unit=8 , f ile= ' /home/noise/sagos/modes/mcoupl . dat ' , 

*  form=' unformatted' , status=' new' ) 
open (unit=14 , f ile= ' /home/noise/sagos/modes/kder . dat ' , 

*  form=' unformatted' , status=' old' ) 

read (4)  w,dx,dy,dz 
rewind  4 

read (14)  ixorig, iyorig 

read (14)  da,na,dr,nr 

write (8)  ixorig, iyorig 

write (8)  da,na,dr,nr 

read (14)  kd 
c 

do  ix=l,nx 

x ( ix) =dx*df loat ( ix- ixorig) 

enddo 
c 

do  iy=l,ny 

y (iy) =dy*df loat (iy- iyorig) 

enddo 
c 

write(6,*)'x  range:  ( ' , x (1) , ' , ' , x (nx) , ' ) ' 

write(6,*)'y  range:  ( ' , y (1) , ' , ' , y (ny) , ' ) ' 

c 

c    computation  of  first  mode  coupling  coefficients 

c 
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c 
c 
c 


do  n=l,nm 
do  m=l,nm 

read (4)  w,dx,dy,dz 
read (4)  nxx,nyy,nzz 
c 

do  ix=l , nx 
do  iy=l,ny 

read(4)  h(ix,iy) 
c 

c    i,  is  the  last  station  in  the  water  column 
c 

i=int (h (ix, iy) /dz)  +  1 
read (4)  cs 
c 

c    calculate  the  total  wavenumber 
c 

do  iz=l,nz 

k (ix, iy, iz) =w/cs (iz) 
enddo 
read (4)  d 
rhol (ix, iy) =d (i) 
rho2 (ix, iy) =d(i+l) 
read (4)  icounter (ix, iy) 
read (4)  kk 
km(ix, iy) =kk (m) 
kn (ix, iy) =kk (n) 
read (4)  efun 
c 

if  (m.eq.n)  then 

ar (m, n, ix, iy) =0 . dO 
aa  (m,  n,  ix,  iy)  =0  .  dO 
else 
c 

c    trapezoid  integration,  to  find  the  integral  part  of 
c    first  coefficient 
c 

sumx>0  .dO 
sumy=0 .dO 

denom=km(ix, iy) -kn (ix, iy) 
c 

do  iz=l,nz 

a=k (ix, iy, iz) *efun (iz, n) * 

*  efun(iz,m) /d(iz) 

sumx=sumx+a*kd (ix, iy, iz, 1) 
sumy = sumy +a*kd (ix, iy, iz, 2 ) 

enddo 
c 

if  (m.eq. 17 .and.n. eq. 16) 

*  write (6 , *) ' ix, iy, sumx, sumy' , ix, iy , 

*  sumx, sumy 
ar (m, n, ix, iy) =4 .dO*sumx*dz/denom 
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c 

c 


c 


aa (m, n, ix, iy) =4 .dO*sumy*dz/denom 
if  (m. eq. 17 . and.n. eq. 16) 

*  write (6, *) ' ix, iy, ar (17 , 16 , ix, iy) ' , 

*  ix, iy, ar (17 , 16 , ix, iy; 
endif 

zbm ( ix, iy) =  efun ( i ,m) 
zbmml (ix, iy) =efun ( i- l,m) 
zbn (ix, iy) =efun ( i, n) 
zbnml (ix, iy) =efun ( i- 1, n) 
enddo 
enddo 

call  part ial ( h , nx , ny , x , y , hpr , hpa , nw , hpx , hpy , c , wk ) 

cal 1  par t  ial ( zbn , nx , ny , x , y , zbnpr , zbnpa , nw , 

*  zbnpx, zbnpy , c, wk) 

if  (m. eq. 17 . and. n. eq. 16)  then 
write (6 , *) ' zbn  =  ' 
write (6 , 10  01) ( ( zbn (ix, iy) , ix=l , nx) , 

*  iy=ny,l,-l) 
write (6 ,*)' zbnpr  =' 

write (6  ,  1001)  ( (zbnpr (ix, iy)  , ix=l, nx) , 

*  iy=ny,l,-l) 
write (6 ,*)' zbnpa  =' 

write (6,1001) ( (zbnpa (ix, iy) , ix=l , nx) , 

*  iy=ny,l,-l) 
1001        format (5 (2x,el2 .4)  ) 

endif 

if  (m.ne.n)  call  sub 1 (rhol , rho2 , dz , zbm, zbn, zbmml , 

*  zbnml , crl, cal, km, kn,m, n, zbnpr, zbnpa, x, y) 

if  (m.eq.n)  call  sub2 (rhol, rho2 , zbm, h, crl, cal , 

*  nx,ny,hpr,hpa,x,y) 

do  ix=l,nx 
do  iy=l,ny 

cr (m,n, ix, iy) =crl (ix, iy) 
ca (m, n, ix, iy) =cal (ix, iy) 
enddo 
enddo 

if  (m. eq. 17 .and.n. eq. 16)  then 

write (6 ,*)' checking  quadrature  :' 

write  (6 , *)  ' ar ( , ix, iy) = ' 

write (6, 100) ( (ar (17, 16 , ix, iy) , ix=l, 5) , iy=l, 5) 

write (6, *)  '  aa  (17,  16 , ix, iy) =' 

write (6, 100)  (  (aa  (17, 16 , ix, iy)  , ix=l, 5)  , iy=l, 5) 
endif 
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c 
c 

c 


do  ix=l,nx 

do  iy=l,ny 

ar (m, n, ix, iy) =ar (m, n, ix, iy) +cr (m, n, ix, iy) 
aa  (m,  n,  ix,  iy)  =aa  (m,  n,  ix,  iy)  +ca  (m,  n,  ix,  iy) 
enddo 
enddo 


rewind  4 
enddo 


do  ix=l,nx 
do  iy=l,ny 

zb (n, ix, iy) =  zbn (ix, iy) 
enddo 
enddo 
c 

enddo 
c 

c    computation  of  second  mode  coupling  coefficients 
c 

do  n=l,nm 
do  m=l,nm 
c 

call  sub3 (rhol, rho2 , zb, h, bl, ar, aa, n,m, hpr, hpa, 

*  zbnpr, zbnpa, x, y) 
c 

do  ix=l,nx 
do  iy=l,ny 

if  (icounter (ix, iy) . It . n. or . 

*  icounter (ix, iy) . It .m)  then 
b  (m,  n,  ix,  iy)  =0  .dO 

else 

b  (m,  n,  ix,  iy)  =bl  (ix,  iy) 
endif 
enddo 
enddo 
c 

enddo 
enddo 
c 

write (6 ,*)' checking  the  mode  coupling  coefficients  :' 
write (6, *) ' ar (16 , 17, ix, iy) =' 

write (6, 100) ( (ar (16 , 17, ix, iy) , ix=l, 5) , iy=l, 5) 
write (6, *) ' aa (16 , 17, ix, iy) = ' 

write (6, 100) ( (aa (16 , 17, ix, iy) , ix=l, 5) , iy=l, 5) 
write (6, *) ' b (16 , 17 , ix, iy) = ' 

write (6, 100) ( (b(16, 17, ix, iy) , ix-1,5) , iy-1, 5) 
100   format (5 (lx,el2 .5) ) 
c 
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write (8)  ar 
write (8)  aa 
write(8)  b 

close (4) 
close (8) 
close (14) 
end 
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************************************************************ 
*  * 

*  This  subroutine  computes  the  vector  mode  coupling  * 

*  coefficient  correction,  due  to  small  bathymetry  * 

*  variations  between  two  different  modes  (m  different  * 

*  than  n.  * 

*  * 
************************************************************ 

c 

subroutine  subl (rhol, rho2 , dz, zbm, zbn, zbmml, zbnml , 

*  crl , cal , km, kn,m, n, zbnpr , zbnpa, x, y) 
c 

parameter  (nx=ll , ny=ll , nw=2*nx*ny+2*max (nx, ny) ) 

implicit  real*8  (a-h,o-z) 

real*8  c (2 , nx, 2 , ny) , wk (nw) , crl (nx, ny) , cal (nx, ny) , 

*  zbm(nx,ny) 
real* 8  zbmml (nx, ny) , zbn (nx, ny) , zbnml (nx, ny) , kn (nx, ny) 
real*8  zbnpz (nx, ny) , zbmpz (nx, ny) , zbnpzpr (nx, ny) , 

*  km  ( nx , ny ) 
real* 8  zbnpzpa (nx, ny) , zbnpzpx (nx, ny) , zbnpzpy (nx, ny) , 

*  zbnpr (nx,ny) 
real* 8  zbnpa (nx, ny) , rhol (nx, ny) , rho2 (nx, ny) 

real*8  x(nx),y(ny) 
c 

do  ix=l,nx 

do  iy=l,ny 

zbnpz (ix, iy) = (zbn (ix, iy) - zbnml (ix, iy) ) /dz 
zbmpz (ix, iy) = (zbm(ix, iy) - zbmml (ix, iy) ) /dz 
enddo 
enddo 
call  partial (zbnpz , nx, ny , x, y , zbnpzpr, zbnpzpa, nw, 

*  zbnpzpx, zbnpzpy, c, wk) 
if  (m. eq. 2 .and.n. eq. 18)  then 

write ( 6 , * ) ' zbnpz (3,3)=' , zbnpz (3,3) 
write ( 6 , * ) ' zbnpzpr (3,3)=' , zbnpzpr (3,3) 
write ( 6 , * ) ' zbnpzpa (3,3)=' , zbnpzpa (3,3) 
write ( 6 , * ) ' zbnpr (3,3)=', zbnpr (3,3) 
write ( 6 , * ) ' zbnpa (3,3)=' , zbnpa (3,3) 
endif 
c 

do  ix=l,nx 
do  iy=l,ny 

r=dsqrt (x(ix) **2  +  y(iy)**2) 

if  (r.lt.l.d-20)  goto  100 

crl (ix, iy) = zbnpzpr (ix, iy) *zbm(ix, iy) * 

*  (I.d0-rho2 (ix, iy) /rhol (ix, iy) ) /rhol (ix, iy)  - 

*  zbmpz (ix, iy) *zbnpr (ix, iy) * 

*  (l.dO/rhol (ix, iy) -I.d0/rho2 (ix, iy) ) 
cal (ix, iy) =zbm(ix, iy) * zbnpzpa (ix, iy) * 

*  (I.d0-rho2 (ix, iy) /rhol (ix, iy) ) / 
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*  (rhol (ix, iy) *r)  -  zbmpz ( ix, iy) * 

*  zbnpa (ix, iy) * 

*  (1. dO /rhol (ix, iy) -I.d0/rho2 (ix, iy) ) /r 

crl  (ix,  iy)  =crl  (ix,  iy)  *2  .6.0/  (kn  (ix,  iy)  -km(  ix,  iy)  ) 
cal  (ix,  iy)  =cal  (ix,  iy)  *2  .6.0/  (kn  (ix,  iy)  -km(ix,  iy)  ) 
100        continue 
enddo 
enddo 

return 
end 
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****************************************************************** 

*  * 

*  This  subroutine  computes  the  vector  mode  coupling         * 

*  coefficient  correction   due  to  small  bathymetry  * 

*  variations,  in  the  case  of  m  equals  n.  * 

*  * 
************************************************************ 
c 

subroutine  sub2 (rhol, rho2 , zbm, h, crl , cal, nx, ny, 

*  hpr,hpa,x,y) 
c 

implicit  real*8  (a-h,o,z) 

real* 8  crl (nx, ny) , cal (nx, ny) , zbm(nx, ny) , rhol (nx, ny) , 

*  rho2  (nx,  ny) 
real*8  hpr(nx,ny) ,hpa(nx,ny) ,x(nx) ,y(ny) 

c 

do  ix=l,nx 
do  iy=l,ny 

r=dsqrt (x(ix) **2  +  y(iy)**2) 
if  (r.lt.l.d-20)  goto  100 
crl(ix,iy)=- (1 .dO/rhol (ix, iy) 

*  -I.d0/rho2 (ix, iy) ) * (zbm(ix, iy) **2) *hpr (ix, iy) 
cal (ix, iy) =- (l.dO/rhol (ix, iy) -I.d0/rho2 (ix, iy) ) * 

*  (zbm(ix, iy) **2) *hpa (ix, iy) /r 
100        continue 

enddo 
enddo 
c 

return 
end 
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*  * 

*  This  subroutine  computes  the  scalar  mode  coupling       * 

*  coefficient  small  bathymetry  changes  included.  * 

*  * 
************************************************************ 


c 


subroutine  sub3 (rhol , rho2 , zb, h, bl , ar , aa, n, m, hpr , hpa, 

*  zbnpr, zbnpa, x, y) 

parameter  (nx=ll , ny=ll , nw=2*nx*ny+2*max (nx, ny) , nm=20) 

implicit  real*8  (a-h,o-z) 

real  *  8  c  ( 2  ,  nx ,  2  ,  ny )  ,  wk  ( nw )  ,  zb  ( run ,  nx ,  ny )  ,  h  ( nx ,  ny )  , 

*  bl (nx, ny) 
real* 8  hpr  (nx,  ny)  ,  hpa  (nx,  ny)  ,  ar  (run,  run,  nx,  ny)  , 

*  aa  ( nm ,  nm ,  nx ,  ny ) 
real*8  zbnpr  (nx,  ny)  ,  zbnpa  (nx,  ny)  ,  ern  (nm,  nx,  ny) 

real* 8  ean (nm, nx, ny) ,sum(nx,ny) ,armn(nx,ny) , 

*  aamn(nx,ny) 
real* 8  arpr(nx,ny) ,arpa(nx,ny) ,arpx(nx,ny) ,arpy(nx,ny) 
real *  8  aapr ( nx , ny )  , aapa ( nx , ny )  , aapx ( nx , ny )  , aapy ( nx , ny ) 
real*8  x(nx) ,y(ny) , erm (nm, nx, ny) , eam(nm, nx, ny) 
real*8  rhol (nx, ny) , rho2 (nx, ny) 

np=l7 
mp=16 

if  (n. eq.np.and.m. eq.mp)  then 

write  (6, *)  ' rhol' 

write (6, 10  0) ( (rhol (ix, iy) , ix-1, 5) , iy=5 , 1,-1) 

write (6, *) ' rho2 ' 

write (6, 10  0) ( (rho2 (ix, iy) , ix=l, 5) , iy=5, 1,-1) 

write (6 , * ) ' hpr' 

write (6, 100) ( (hpr (ix, iy) , ix=l, 5) , iy=5, 1,-1) 

write(6, * ) ' zb ( ' ,m, ' , . . . ' 

write (6, 100)  ( (zb (m, ix, iy)  , ix=l, 5)  , iy=5 ,1, - 1) 

write(6,*) ' zb ( ' , n, ' , . . . ' 

write (6, 100) ( (zb(n, ix, iy) , ix=l, 5) , iy=5, 1,-1) 
endif 

do  ix=l,nx 

do  iy=l,ny 

sum(ix, iy) =0 . dO 

enddo 
enddo 

do  1=1 , nm 

do  ix=l,nx 
do  iy=l,ny 

r=dsqrt (x(ix) **2  +  y(iy)**2) 
if  (r.lt.l.d-20)  goto  110 
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ern (1 , ix, iy) = . 5d0*ar (n, 1 , ix, iy)  +  hpr(ix,iy)* 

*  (l.dO/rhol (ix, iy) -I.d0/rho2 (ix, iy) ) * 

*  zb (n,  ix, iy) *zb (1 , ix, iy) 

ean (1 , ix, iy) = . 5d0*aa (n, 1 , ix, iy)  +  hpa(ix,iy)* 

*  (l.dO/rhol (ix, iy) -I.d0/rho2 (ix, iy) ) * 

*  zb (n, ix, iy) *zb (1 , ix, iy) /r 

errnd , ix, iy) = . 5d0*ar (m, 1 , ix, iy)     +   hpr(ix,iy)* 

*  (l.dO/rhol (ix, iy) -I.d0/rho2 (ix, iy) ) * 

*  zb (m, ix, iy) *zb (1 , ix, iy) 

earn (1 , ix, iy) = . 5d0*aa (m, 1 , ix, iy)     +   hpa(ix,iy)* 

*  (l.dO/rhol (ix, iy) -I.d0/rho2 (ix, iy) ) * 

*  zb (m, ix, iy) *zb (1 , ix, iy) /r 
110        continue 

enddo 
enddo 
enddo 

i 

do  ix=l,nx 
do  iy=l,ny 
do  1=1, nm 

sum(ix, iy) =sum(ix, iy)  + 

*  ern (1 , ix, iy) *erm(l , ix, iy) + 

*  ean (1 , ix, iy) *eam(l , ix, iy) 
enddo 

enddo 
enddo 

p 

if  (n. eq.np.and.m. eq.mp)  then 

write (6 , *) ' sum  for  im, in= ' ,mp, np 
write(6,100) ( (sum(ix, iy) , ix=l, 5) , iy=5, 1,-1) 
100  format (5 (lx,el2.3) ) 

endif 

do  ix=l,nx 
do  iy=l,ny 

armn  (ix,  iy)  =ar  (m,  n,  ix,  iy) 
aamn  (ix,  iy)  =aa  (m,  n,  ix,  iy) 
enddo 
enddo 

i 

call  partial (armn, nx, ny, x, y, arpr, arpa, nw, 

*  arpx, arpy, c, wk) 
call  partial (aamn, nx, ny,x, y, aapr, aapa, nw, 

*  aapx, aapy, c, wk) 

t 

do  ix=l,nx 
do  iy=l,ny 

r=dsqrt (x(ix) **2  +  y(iy)**2) 

if  (r.lt .l.d-20)  then 

bl ( ix, iy) = . 5d0*arpr (ix, iy)  -  sum(ix,iy)  - 

*  (l.dO/rhol (ix, iy) -I.d0/rho2 (ix, iy) ) *zb (m, ix, iy) * 
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*  (hpr (ix, iy) *zbnpr (ix, iy)  ) 
else 

bl (ix, iy) = . 5d0*arpr (ix, iy) + . 5d0*ar (m, n, ix, iy) /r+ 

*  . 5*d0*aapa ( ix, iy) /r  -  sum(ix,iy)  - 

*  (l.dO/rhol (ix, iy) - 

*  I.d0/rho2 (ix, iy) ) *zb (m, ix, iy) * 

*  (hpr (ix, iy) *zbnpr (ix, iy) + 

*  hpa ( ix, iy) *zbnpa (ix, iy) /r**2 ) 
endif 

enddo 
enddo 

return 
end 
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************************************************************ 

*  * 

*  This  subroutine  computes  the  partial  derivatives  with     * 

*  respect  to  range  and  azimuthal  angle  of  a  given  function  * 

*  f (x,y) .  It  uses  a  bicubic  spline  to  calculate  the         * 

*  cartesian  derivatives  and  then  perform  a  coordinate      * 

*  transformation  to  cylindrical  coordinates.  * 

*  * 
*********************************************************** 

c 

subroutine  partial (f , nx, ny,x,y, fpr, fpa, nw, fpx, fpy, 

*  c , wk ) 

c 

implicit  real*8  (a-h,o-z) 

real*8  f (nx, ny) , x (nx) , y (ny) , fpr (nx, ny) , fpa (nx, ny) 

real* 8  wk (nw) , fpx (nx, ny) , fpy (nx, ny) , c (2 , nx, 2 , ny) 


c 

c 

c 
c 


external  ibcccu 

ic=nx 

pi=dacos ( - 1 . dO) 

call  ibcccu (f , x, nx,y, ny, c, ic, wk, ier) 

do  ix=l,nx 

do  iy=l,ny 

fpx(ix, iy) =c (2, ix, 1, iy) 
fpy (ix, iy) =c (1, ix,2, iy) 
if  (x(ix) .eq.O.dO)  then 

theta=dsign(y (iy) , l.dO) *pi/2 .dO 
else 

theta=datan (y ( iy) /x ( ix) ) 
endif 

r=dsqrt (x(ix) **2+y (iy) **2) 
fpr (ix, iy) =fpx (ix, iy) *dcos (theta) + 

*  fpy (ix, iy) *dsin(theta) 

fpa (ix, iy) =2  fpx (ix, iy) *r*dsin (theta)  + 

*  fpy (ix, iy) *r*dcos (theta) 
enddo 

enddo 

return 
end 
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